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Abstract 

The Lambert W(x) function and its possible applications in physics are presented. The actual numerical implementa- 
tion in C++ consists of Halley's and Fritsch's iterations with initial approximations based on branch-point expansion, 
^> asymptotic series, rational fits, and continued-logarithm recursion. 
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^ Program summary 

i 1 Program title: LambertW 

1/1 Catalogue identifier: AENC_vl_0 

Program summary URL: http://cpc.cs.qub.ac.uk/summaries/AENC_vl_0.html 
• Program obtainable from: CPC Program Library, Queen's University, Belfast, N. Ireland 
q Licensing provisions: GNU General Public License version 3 
'— 'No. of lines in distributed program, including test data, etc.: 1335 

| No. of bytes in distributed program, including test data, etc.: 25283 
^ Distribution format: tar.gz 

Programming language: C++ (with suitable wrappers it can be called from C, Fortran etc.), the supplied command-line 
C*~) utility is suitable for other scripting languages like sh, csh, awk, perl etc. 
t • Computer: All systems with a C++ compiler. 

Operating system: All Unix flavors, Windows. It might work with others. 
RAM: Small memory footprint, less than 1 MB 
O Classification: 1.1, 4.7, 11.3, 11.9. 

Nature of problem: Find fast and accurate numerical implementation for the Lambert W function. 
Solution method: Halley's and Fritsch's iterations with initial approximations based on branch-point expansion, 
asymptotic series, rational fits, and continued logarithm recursion. 
^ Additional comments: Distribution file contains the command-line utility lambert-w. Doxygen comments, included in 
!-h the source files. Makefile. 

Running time: The tests provided take only a few seconds to run. 
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1. Introduction 

The Lambert W function is defined as the inverse 
function of the 

i4i e x (1) 



mapping and thus solves the 



(2) 



equation. This solution is given in the form of the Lam- 
bert W function, 

y = w(x), (3) 

and according to Eq. (2) it satisfies the following relation, 

W(x) e w W = x. (4) 

Since the mapping in Eq. (1) is not injective, no unique 
inverse of the xe x function exists (except at the mini- 
mum). As can be seen in Fig. 1, the Lambert W function 
has two real branches with a branching point located at 
( — e , —1). The bottom branch, W_i(x), is defined in 
the interval x E \—e , 0] and has a negative singular- 
ity for x -4- _ . The upper branch Wq(i) is defined for 
x 6 [-e _1 , °°]. 

The earliest mention of the problem of Eq. (2) is at- 
tributed to Euler [1]. However, Euler himself credited 
Lambert for his previous work in this subject, Lambert's 
transcendental equation [2]. The W(x) function started 
to be named after Lambert only recently, in the last 20 
years or so; nevertheless, the first usage of the letter W 
appeared much earlier [3]. 

Recently, the W(x) function amassed quite a follow- 
ing in the mathematical community [4] . Its most faithful 
proponents are suggesting elevating it among the present 
set of elementary functions, such as sin(x), cos(x), exp(x), 
ln(x), etc. The main argument for doing so is the fact 
that W is the root of the simplest exponential polynomial 
function given in Eq. (2). We will acknowledge these ef- 
forts by strict usage of a roman symbol W as its name. 

While the Lambert W function is called LambertW in 
the mathematics software tool Maple [5] and lambertw in 
the Matlab programming environment [6], in the Math- 
ematica computer algebra framework this function [7] is 
implemented under the name ProductLog (in the recent 
versions an alias LambertW is also supported). In open 
source format the Lambert W function is available in 
the special-functions part of the GNU Scientific Library 
(GSL) [8], unfortunately implemented using only the slower 
Halley's iteration (see discussion in Section 6). 

There are numerous, well documented applications 
of W(x), certainly abundant in mathematics (like lin- 
ear delay-differential equations [9]), numerics [10], com- 
puter science [11] and engineering [12], but surprisingly 
many also in physics, just to mention a few (without 
preference): quantum mechanics (solutions for double- 
well Dirac-delta potentials [13]), quantum statistics [14], 
general relativity (solutions to (l+l)-gravity problem [15], 
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Figure 1: The two branches of the Lambert W function, W_j (x) in blue 
and Wo(i') in red. The branching point at (— e , — 1) is indicated with 
a green dash. 



inverse of Eddington-Finkelstein (tortoise) coordinates [16]), 
statistical mechanics [17], fluid dynamics [18], optics [19] 
etc. 

The main motivation for the implementation in this 
work comes from the research in cosmic ray physics where 
it has been in use already for several years [20] and new 
interesting applications are appearing frequently [21]. In 
the next sections let us describe two new examples that 
arise from this field of astrophysics. 

1.1. Inverse of the May al function 

The Moyal function and the related distribution was 
proposed as a good approximation for the more com- 
plicated Landau distribution [22] used in descriptions of 
energy loss of charged particles passing through matter 
[23]. The un-normalized Moyal function is defined as 



M(x) = exp 



■Ux + e~ x ) 



(5) 



and can be seen in Fig. 2 (top). 

Its inverse can be written in terms of the two branches 
of the Lambert W function, 



M^}(x) = W _i(-x z ) -21nx. 



(6) 



In experimental astrophysics the Moyal function is 
used for phenomenological recovery of the saturated sig- 
nals from photomultipliers [24], where the largest parts 
of the saturated signals are obscured by the limited range 
of the analog-to-digital converters. 

1.2. Inverse of the Gaisser-Hillas function 

In astrophysics the Gaisser-Hillas function is used to 
model the longitudinal particle density in a cosmic-ray 
air shower [25]. We can show that the inverse of the 
three-parametric Gaisser-Hillas function, 



G(X; Xq, X max , A) 



X-Xp 
-X 



X n 



A 



exp 



X n 



X s 



(7) 
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Figure 2: Top: The Moyal function M(x). Bottom: A family of one- 
parametric Gaisser-Hillas functions g(x; x max ) for x max in the range 
from 1 to 10 with step 1. 



is intimately related to the Lambert W function. Using 
rescale substitutions, 



x-x 

A 



and 



X n 



A 



(8) 



the Gaisser-Hillas function is modified into a function of 
one parameter only, 



exp(x n 



x). 



(9) 



The family of one-parametric Gaisser-Hillas functions is 
shown in Fig. 2 (bottom). The problem of finding an 
inverse, 

g(x; x max ) = y (10) 
for < y 1, can be rewritten into 



exp 



-y 



1 / %ax ^ 1 



(11) 



According to the definition (2), the two (real) solutions 
for x are obtained from the two branches of the Lambert 
W function, 



Xl,2 



XmaxWo-U-y 



1/Xn 



Note that the branch —1 or simply chooses the right or 
left side relative to the maximum, respectively. 

The Gaisser-Hillas function is intimately related to 
the gamma distribution which was successfully used some- 
what earlier [26] in an approximate description of the 
mean longitudinal profile of the energy deposition in 
electromagnetic cascades. It is trivial to show that the 
inverses of the gamma and inverse-gamma distributions 
[27] are expressible in terms of the Lambert W function 
as well. 



2. Numerical methods 

Before describing the actual implementation let us re- 
view some of the possible numerical and analytical ap- 
proaches to find W(x). 

2.1. Recursion 

For x > and W(x) > we can take the natural 
logarithm of Eq. (4) and rearrange it into 



W(x) = lnx-lnW(x). 



(13) 



From here it is clear that an analytical expression for 
W(x) will exhibit a certain degree of self similarity. The 
W(x) function has multiple branches in the complex do- 
main. Due to the x > and W(x) > conditions, 
Eq. (13) represents the positive part of the Wrj(x) prin- 
cipal branch and in this form it is suitable for evaluation 
when Wo(x) > 1, i.e. when x > e. 

Unrolling the self-similarity in Eq. (13) as a recursive 
relation, the following curious expression for Wq(x) is 
obtained, 



Wo(x) = lnx — ln(lnx — ln(lnx — ...)), 
or in the shorthand of a continued logarithm, 

W (x)=ln : 



In 



(14) 



(15) 



In 



The above expression clearly has the form of successive 
approximations, the final result given by the limit, if it 
exists. 

For x < and W(x) < we can multiply both sides 
of Eq. (4) with —1, take logarithm, and rewrite it into a 
similar expansion for the W_x(x) branch, 



(12) 



W(x) =ln(-x) -ln(-W(x)). 
Again, this leads to the similar recursion, 

W_ a (x) = ln(-x) -ln(-(ln(-x) - ln(- 
or as a continued logarithm, 

W_ x (x) = In- 



to. - 



(16) 



•)))/ (17) 



(18) 



■ln- 



3 



For this continued logarithm we will use the symbol R_\ (x) 
where n denotes the depth of the recursion in Eq. (18). 
Starting from yet another rearrangement of Eq. (4), 



The relative difference can be expressed as 



W(x) 



exp W(x)' 



we can obtain a recursion relation for the — e 1 < x < e 
part of the principal branch Wq(x) < 1, 



W (at) 



exp 



(20) 



exp — 



2.2. Halley's iteration 

We can apply Halley's root-finding method [28] to de- 
rive an iteration scheme for f(y) = W(y) — X by writing 
the second-order Taylor series 

f(y) = f{yn) +f{yn) (y - y«) + \f"{y n ) (y - y«) 2 + ■ ■ ■ 

(21) 

Since root y of f(y) satisfies f(y) — we can approxi- 
mate the left-hand side of Eq. (21) with and replace y 
with y„+i. Rewriting the obtained result into 



y»i+i — Vn 



f(lfn) 



f'(yn) + \f"{y n ) (y n+ i-y n ) 



(22) 



and using Newton's method y n +\ ~yn = —f(yn)/f"(yn) 
on the last bracket, we arrive at the expression for Hal- 
ley's iteration for the Lambert W function 



Vf n 



tn 



where 



t„ = W n e w » -x, 



W„+2 



2(W» + 1)' 
u„ = (W„ + l)e w «. 



(23) 

(24) 
(25) 
(26) 



This method is of the third order, i.e. having W n = 
W(x) + 0(e) will produce W„ +1 = W(x) + 0{e 3 ). Sup- 
plying this iteration with a sufficiently accurate first ap- 
proximation of the order of O(10 -4 ) will thus give a 
machine-size floating point precision (9(10 -16 ) in at least 
two iterations. 

2.3. Fritsch's iteration 

For both branches of the Lambert W function a more 
efficient iteration scheme exists [29], 



W n+ i = W„ (! + £„), 



(27) 



where e n is the relative difference of successive approxi- 
mations at iteration n, 



Wn+l ~ W n 



(28) 



£n = 



1 + W„J \q„ -2z„ 



(19) where 



z„ = In — - W„, 
W„ 

q n =2(l + W n )(l + W n + lz n) 



(29) 

(30) 
(31) 



The error term in this iteration is of a fourth order, i.e. 
with W n = W(x) + 0(e n ) we obtain W, I+1 = W(x) + 
(9(4). 

Supplying this iteration with a sufficiently reason- 
able first guess, accurate to the order of (9(10~ 4 ), will 
therefore deliver machine-size floating point precision 
O(10~ 16 ) in only one iteration and excessive (9(10 -64 ) 
in two! The main goal now is to find reliable first order 
approximations that can be fed into Fritsch's iteration. 
Due to the lively landscape of the Lambert W function 
properties, the approximations will have to be found in 
all the particular ranges of the function's behavior. 

3. Initial approximations 

The following section deals with finding the appro- 
priate initial approximations in the whole definition range 
of the two branches of the Lambert W function. 

3.1. Branch-point expansion 

The inverse of the Lambert W function, W _1 (y) = 
yeV ', has two extrema located at W~ 1 ( — 1) = — e -1 and 
W _1 (— oo ) = 0~. Expanding W _1 (y) to the second order 
around the minimum at y = —1 we obtain 



W-\y) 



(y + i) ; 

2e 



(32) 



The inverse W _1 (y) is thus approximated in the lowest 
order by a parabolic term implying that the Lambert W 
function will have square-root behavior in the vicinity of 
the branch point x = — e _1 , 



W_i, (*) 



-1T\/2(1 



ex 



(33) 



To obtain the additional terms in Eq. (33) we proceed 
by defining an inverse function, centered and rescaled 
around the minimum, i.e. f(y) — 2(eW _1 (y — 1) + 1). 
Due to this centering and rescaling, the Taylor series around 
y = becomes particularly neat, 



/(y) = y 2 + iy 3 + iy 4 + ^y 5 



(34) 



Using this expansion we can derive coefficients [30] of 
the corresponding inverse function 



f- l {z) = 1+W 

- vV2 _ 



2e 



11-3/2 _ 43__2 



72" 



540" 



(35) 
(36) 
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Figure 3: Successive orders of the branch-point expansion for W_i (x) 
on the top and Wo(x) on the bottom. 

From Eq. (35) we see that z = 2(1 + ex). Using p±(x) = 
±-\/2(l + ex) as an independent variable we can write 
this series expansion as 

Vf-ljo(x) « BL"l /0 (x) = f>p^(x), (37) 

!=0 

where the lowest few coefficients h are 
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and more of them are given in the accompanying code 
(see Fig. 3 for succesive orders of the series). 

3.2. Asymptotic series 

Another useful tool is the asymptotic expansion [31]. 
Using 



A(a,b) = a - b + J^J^C km a 

k m 



-k—m—l r.m+1 



(38) 



with Q- m related to the Stirling numbers of the first kind, 
the Lambert W function can be expressed as 



Vf-uo(x) = A(]n(Tx), ln(Tln(Tz))), 



(39) 



with a = Inx, b = lnlnx for the Wq branch and a = 
ln(— x), b = ln(— ln(— x)) for the W_i branch. The first 
few terms are 



A(a,b) = a-b^ h 



b b(-2 + b) b(6-9b + 2b 2 ) 



2a 2 



6« 3 



+ 



b(-12 + 36b-22b 2 + 3b 3 ) 

12a 4 ~ + 

b(60 - 300b + 350b 2 - 125b 3 + 12b A ) 
60o5 



(40) 



3.3. Rational fits 

A useful quick-and-dirty approach to the functional 
approximation is to generate a large enough sample of 
data points {iVj e w ', w{\, which evidently lie on the Lam- 
bert W function. Within some appropriately chosen range 
of Wj values, the points are fitted with a rational approx- 
imation 

Yd ajX 1 



Q(x) 



(41) 



Ydkx 1 ' 

varying the order of the polynomials in the nominator 
and denominator, and choosing the one that has the low- 
est maximal absolute residual in a particular interval of 
interest. 

For the Wo(x) branch, the first set of equally-spaced 
Wj components was chosen in a range that produced 
u>i e w ' values in an interval [—0.3, 0] . The optimal rational 
fit turned out to be 



Qo(x) 



1 + a\x + aj_x 2 + a$x 3 + a^x 
1 + b x x + b 2 x 2 + b 3 x 3 + b^x 4 



(42) 



where the coefficients 1 for this first approximation Qfo (x) 
are 
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For the second fit of the Wo(x) branch a if, range 
was chosen so that the if, e Wi values were produced in 
the interval [0.3, 2e\, giving rise to the second optimal 

f2l 

rational fit Q 1 (x) of the same form as in Eq. (42) but 
with coefficients 
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For the W_i(x) branch a single rational approxima- 
tion of the form 



Q-i(x) 



«0 + a\x + aix 1 



1 + b\x + b2X 2 + b^x 3 + + b^x 5 
with the coefficients 



(43) 



1 The ' symbol in coefficient values denotes truncation in this pre- 
sentation; the full machine-size sets of decimal places are given in the 
accompanying code. 
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is enough. 

4. Implementation 

To quantify the accuracy of a particular approxima- 
tion W(x) of the Lambert function W(x) we can intro- 
duce a quantity A(x) defined as 

A(x) = log 10 | W(x) | - log 10 | W(x) - W(x)|, (44) 

so that it gives us a number of correct decimal places the 
approximation W(x) is producing for some value of the 
parameter x. 

In Fig. 4 all approximations for the Wo(x) mentioned 
above are shown in the linear interval [— e _1 , 0.3] on the 
left and the logarithmic interval [0.3, 10 5 ] on the right. 
For each of the approximations an optimal interval is 
chosen so that the number of accurate decimal places 
is maximized over the whole definition range. For the 
Wo(x) branch the resulting piecewise approximation 

191 



Wo(ac) 



TO 



-e 1 < x < a 



; a ^ x < b 
; b ^ x < c 
; c ^ x < oo 



(45) 



with a = -0.323581', b = 0.145469', and c = 8.706658', is 
accurate in the definition range [— e , 7] to at least 5 dec- 
imal places and to at least 3 decimal places in the whole 

definition range [— e _1 , oo]. Note that Bq (x) comes from 

Eq. (37), Q$ (x) and (x) are from Eq. (42), and A (x) 
is from Eq. (40). 

The accuracy of the final piecewise approximation 
Wo(i) is shown in Fig. 5 (black line). Using this ap- 
proximation, a single step of Halley's iteration (red) and 
Fritsch's iteration (blue) is performed and the resulting 
number of accurate decimal places is shown. As can 
be seen from the plots, both iterations produce machine- 
size accurate floating point numbers in the whole defini- 
tion interval except for the interval between 6.5 and 190 
where the Halley's method requires another step of the 
iteration. For that reason we have decided to use only 
(one step of) Fritsch's iteration in the C++ implementa- 
tion of the Lambert W function. 

In Fig. 6 (top) the same procedure is shown for the 
W_i(x) branch. The final approximation 



W_i(x) 




^ x < a 



; a ^ x < b 
;b^x<0 



(46) 




20 



15 



5^- 



10 



100 



1000 



10 4 



10 3 



Figure 4: Combining different approximations of Wo(i) into a final 
piecewise function. The number of accurate decimal places A(x) is 
shown for two ranges, linear interval [— e -1 , 0.3] (top) and logarith- 
mic interval [0.3, 10 5 ] (bottom). The approximations are branch-point 
expansion B^\x) from Eq. (37) (blue), rational fits Q^\x) and q|, 2 '(x) 
from Eq. (42) in black and red, respectively, and asymptotic series Ag(x) 
from Eq. (40) (green). 



with a = -0.302985' and b = -0.051012', is accurate 
to at least 5 decimal places in the whole definition range 

[-e _1 , 0]. Note that (x) is taken from Eq. (37), Q_j (x) 
is from Eq. (43), and R^^x) is from Eq. (18). 

In Fig. 6 (bottom) the combined approximation W_j (x) 
is shown (black line). The values after one step of Hal- 
ley's iteration are shown in red and after one step of 
Fritsch's iteration in blue. Similarly as for the previous 
branch, Fritsch's iteration turns out to be superior, yield- 
ing machine-size accurate results in the whole definition 
range, while Halley's iteration is accurate to at least 13 
decimal places. 

5. Source availability, installation and usage 

The most recent version of the sources of this imple- 
mentation with some additional material and examples 
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Figure 5: Final values of the combined approximation Wq(x) (black) 
from Fig. 4 after one step of Halley's iteration (red) and one step of 
Fritsch's iteration (blue). 

are available at http://www.ung.si/~darko/LambertW/ 
and are released under the GPL license. 

Apart from the special functions in GSL [8], this is 
the only freely available implementation of the Lambert 
W function in C++ and to the best of our knowledge the 
only implementation using the superior Fritsch's version 
of the iteration. 

The supplied C++ code implements the following set 
of functions 2 

• double LambertWApproximation<£» (double x) ; 

• double LambertW<£>> (double x) ; 

• double LambertW(int branch, double x) ; 

where b in the first two functions should be replaced with 
the compile-time branch integer value, e.g. LambertW<-l> (x 
or LambertW<0> (x) . Apart from the slightly increased ef- 
ficiency the main reason for implementing the first two 
functions with the branch b as a compile time parameter 



"Which can be found in the files LambertW . h and LambertW . cc. 
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Figure 6: Top: Approximations of the W_i(x) branch. The branch- 
point expansion B^\(x) is shown in blue, the rational approximation 

Q_l (x) in black, and the logarithmic recursion in red. Bottom: The 
combined approximation is accurate to at least 5 decimal places in the 
whole definition range. The results after applying one step of Halley's 
iteration are shown in red and after one step of Fritsch's iteration in 
blue. 

is to force the users to think about the two possible so- 
lutions to the problem in Eq. (2). Just as for solutions to 
the quadratic equation where the ± sign has to be chosen 
based on the desired solution, the Lambert W function 
offers two possibilities that need careful consideration. 

The initial approximations Wj,(x) from Eqs. (46) and 
(45), that are used to kick-start the iterations, are also 
directly available as LambertWApproximation<£>>(x), as 
they might be useful in applications for which it is suffi- 
cient to have a limited number of accurate decimal places 
(see the discussion above). 

The supplied code does not need any kind of spe- 
cial installation procedures. In the case that your analy- 
sis needs an evaluation of the Lambert W function, the 
two source files, LambertW. h and LambertW. cc, should 
be simply bundled with your project structure and com- 
piled with a suitable C++ compiler. 

The source distribution also includes a command-line 
utility implemented by the lambert-w . cc source file. The 
included Makefile can, with the request make lambert-w, 
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Figure 7: Execution speedup fLj/f/ relative to the GSL implementa- 
tion [8] for the Wg(x) branch (red) and the W_i(x) branch (blue). For 
x > 2 the ratio slowly decreases and is ~ 2 for x > 8. Time of the 
respective driver loops and function calls was subtracted in order to 
measure only differences between the pure numerical parts of the two 
implementations (see text for details). 

build an executable command. It can be invoked through 
a shell as . /lambert-w [branch] x, taking an optional 
branch, number (0 by default) and a parameter x. The 
output of the command is equivalent to the Wb r anch( x ) 
return value and can thus be simply used in shell scripts 
(sh, bash, or csh) or other programming languages with 
easy access to shell invocations (awk, perl etc.). 

Also included in the distribution is a test suite which 
can perform a correctness check on your build by com- 
paring obtained and expected return values of the Lam- 
bert W function on your system. It is invoked by the com- 
mand make tests. Any potential discrepancies larger 
than the double machine precision (> 10~ 14 ) will be re- 
ported in the output. 

6. Timing 

We have decided to compare the execution time of 
our code relative to the GNU GSL library (implemented 
in the C language) since comparisons to interpreted code 
(like Maple, Matlab or Mathematica) would, besides com- 
mon availability problems, not be fair in terms of speed. 

In Fig. 7 relative speedup, f gs i/ i, is shown as a func- 
tion of the Lambert W parameter x. Timing accuracy 
of several % was achieved by running 3 000 000 function 
calls in a loop, taking special care that the compiler did 
not optimize code away by slightly modifying x on each 
call and then gathering all results into a summed variable 
reported at the end. For each of the two implementations 
an identical code was also run with the Lambert W func- 
tion call replaced with a simple identity function (just 
returning its input parameter) in order to estimate the 
overhead of the surrounding timing code. This f verhead is 
then consequently subtracted from the time of the Lam- 
bert W runs t, giving an approximation to the time taken 



by the pure function call, f' = i — f over head- The ratio 
fg gl / 1' is then plotted in Fig. 7. 

As can be clearly seen from Fig. 7, our implementa- 
tion is at least 2x faster than GSL for a broad range of 
input parameters x, but the largest efficiency gains (up 
to 5 x ) are observed in the ranges where rational fits Qj, 
from Eqs. (45) and (46) are used. Although Fritsch's iter- 
ation is in general more complex than Halley's, at the end 
it pays off, yielding machine-size accuracy with a single 
step where Halley's might need more, also due to poor 
initial approximations used in GSL. GSL performs better 
only in the small regions where branch-point and asymp- 
totic expansions are used without the consequent Hal- 
ley's iteration refinements. The comparisons were made 
on the Ubuntu 12.04 x86_64 Linux operating system run- 
ning on a 2.2 GHz AMD Opteron 275 processor, using 
the GCC 4.6.3 compiler with optimization option -02. 

7. Conclusions 

We have shown that Fritsch's iteration scheme cou- 
pled with accurate initial approximations can deliver sig- 
nificant efficiency gains in the numerical evaluation of 
the real branches of the Lambert W function. The open- 
source release of our C++ implementation is suitable for 
inclusion in various analysis packages used in all fields 
of physics. 
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